// RosenbrockFunction.cpp: implementation of the RosenbrockFunction class.
// when run on 20 dimensions, the output in zakharov1.out is interesting
//////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <fstream>
#include <string>

#include "rosenbrock_function.h"
#include "DataHolder.h"



//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

RosenbrockFunction::RosenbrockFunction(int i_NInput)
{
  NInput = i_NInput;
  nCall = 0;

  // Create 3 arrays
  newA(mDerivA, NInput);
  newA(dDerivative, NInput);  
  newA(mDeriv, NInput);
  newA(mDerivative, NInput);

  Result.Size(NInput);
}


RosenbrockFunction::~RosenbrockFunction()
{
  deleteA(mDerivA);
  deleteA(dDerivative);
  deleteA(mDeriv);
  deleteA(mDerivative);
}


// takes a point and if it is beyong boundary move the point
// into the boundary box
void RosenbrockFunction::checkBoundary(vector<double> &Input)
{
  for (unsigned int i = 0; i < Input.size(); i++)
  {
    // reflective boundary
    /*while( (Input[i] >= 10.0) || (Input[i] <= -5.0) )
    {
      if ( Input[i] >= 10.0 )
      {
        Input[i] = 10.0 - fabs(Input[i]-10.0);
      }
      else if ( Input[i] <= -5.0 )
      {
        Input[i] = -5.0 + fabs(-5.0-Input[i]);
      }
    }*/

    // periodic boundary
    if ( Input[i] >= 10.0 )
    {
      Input[i] = -5.0 + fmod(Input[i],10.0);
    }
    else if ( Input[i] <= -5.0 )
    {
      Input[i] = 10.0 + fmod(Input[i],5.0);
    }

    // hard boundaries
    /*if ( Input[i] >= 10.0 )
    {
      Input[i] = 10.0;
    }
    else if ( Input[i] <= -5.0 )
    {
      Input[i] = -5.0;
    }*/
  }
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void RosenbrockFunction::checkBoundary(double *Input)
{
  for (int i = 1; i <= NInput; i++)
  {
    // periodic boundary
    if ( Input[i] >= 10.0 )
    {
      Input[i] = -5.0 + fmod(Input[i],10.0);
    }
    else if ( Input[i] <= -5.0 )
    {
      Input[i] = 10.0 + fmod(Input[i],5.0);
    }
  }
}


//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

double RosenbrockFunction::Value(double *Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 0.0;

  // Makes answer equal to it self added to
  // the sum of x^2/4000 for each dimension
  for (int m = 1; m < NInput; m++)
  {
    answer += ((Input[m]*Input[m])-Input[m+1])*((Input[m]*Input[m])-Input[m+1]);
  }

  answer *= 100;

  // Makes danswer equal to itself multiplied
  for(int n=1; n<NInput; n++)
  {
    danswer += (Input[n]-1)*(Input[n]-1);
  }
  
  answer += danswer;

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// Calculates/builds the function value
double RosenbrockFunction::Value(vector<double> &Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);


  // Initialise vars answer and danswer (both local)
  double answer = 0.0;
  double danswer = 0.0;

  // Constructing the equation
  for (int m = 0; m < (NInput-1); m++)
  {
    answer += ((Input[m]*Input[m])-Input[m+1])*((Input[m]*Input[m])-Input[m+1]);
  }

  answer *= 100;

  for(int n=0; n<(NInput-1); n++)
  {
    danswer += (Input[n]-1)*(Input[n]-1);
  }
  
  answer += danswer;

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// Calculates/builds the derivative (used in CGM and HMC)
void RosenbrockFunction::BuildDerivative(double *Input, double *Derivative)
{
  // NInput = dimension of function

  // Constructing the equation
  for(int r=1; r<=NInput; r++)
  {
    Derivative[r] = 0.0;
  }


  for(int i=1; i<=(NInput-1); i++)
  {
    mDeriv[i] = 2*(Input[i]-1);
    dDerivative[i] = (400*Input[i])*((Input[i]*Input[i])-Input[i+1]);
  }
  dDerivative[NInput] = 0.0;
  mDeriv[NInput] = 0.0;

  mDerivA[1] = 0.0;
  for(int m=2; m<=NInput; m++)
  {
    mDerivA[m] = (-200)*((Input[m-1]*Input[m-1])-Input[m]);
  }

  for(int n=1; n<=NInput; n++)
  {
    Derivative[n] = mDerivA[n] + dDerivative[n] + mDeriv[n];
  }
  
  // stores the derivative as an array
  // first element = derivative with respect to the first var
  for(int j=1; j<=NInput; j++)
  {
    mDerivative[j] = Derivative[j];
  }
}

// Used to test the function is working
void RosenbrockFunction::test_function()
{
  ofstream out("output/test_of_rosenbrock.txt");
  
  out << "NUMBER_OF_DIMENSIONS: " << NInput << endl 
    << "      Subscript[x,k]      |     Value    |         Derivative[k]" << endl << endl;

  int NRun = 3;
  double *mInput;
  double *mDerivArray;
  newA(mInput, NInput);
  newA(mDerivArray, NInput);

  for(int run=1; run<=NRun; run++)
  {
    for(int n=1; n<=NInput; n++)
    {
      mInput[n]=(15.0*ran1())-5.0;
    }

    for(int k=1; k<=NInput; k++)
    {
      out << setw(11) << mInput[k];
    }
    out << setw(16) << Value(mInput);

    BuildDerivative(mInput, mDerivArray);

    for(int l=1; l<=NInput; l++)
    {
      out << setw(16) << mDerivative[l];
    }

    out << "\n\n";
    
  }

  deleteA(mInput);
  deleteA(mDerivArray);
  out.close();
}

// used to print the solution (this is the old version, the new version is below)
void RosenbrockFunction::print_sol(double *solution, double funcVal, char filename[])
{
  ofstream out(filename);	

  out << "Number of dimensions: " << NInput << endl;
  out << "Final Rosenbrock value = " << funcVal << endl;
  out << "Solution coordinates:" << endl;

  for (int i = 1; i <= NInput; i++)
  {
    out << solution[i] << "  ";
  }
 
  out << endl;

  out.close();
}

// used to print the solution to a file
void RosenbrockFunction::print_sol(char filename[])
{
  ofstream out(filename);	

  // getBestValue() returns the best value stored in the tracker
  out << "Number of dimensions: " << NInput << endl;
  out << "Final Rosenbrock value = " << Result.getBestValue() << endl;
  out << "Solution coordinates:" << endl;
  
  // prints out the solution coordinates to a file
  Result.printSolution(out);

  out.close();
}

// distance from know global minimum which is located at the origin
double RosenbrockFunction::distanceFromGlobalMinimum(double *solution)
{
  double toReturn = 0.0;

  for (int i = 1; i <= NInput; i++)
  {
    toReturn += (solution[i]-1.0)*(solution[i]-1.0);
  }

  double toReturnS = sqrt(toReturn);
  return toReturnS;
}